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The canonical partition function approach was designed to avoid the overlap problem that affects 
the lattice simulations of nuclear matter at high density. The method employs the projections 
of the quark determinant on a fix quark number sector. When the quark number is large, the 
evaluation of the projected determinant becomes numerically unstable. In this paper a different 
evaluation method based on expanding the determinant in terms of loops winding around the 
lattice is studied. We show that this method is stable and significantly faster than our original 
algorithm. This greatly expands the range of quark numbers that we can simulate effectively. 
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1. Introduction 

To simulate QCD at finite density, the usual approach based on the grand canonical partition 
function is to split the fermionic determinant into two parts: a real and positive part used to gen- 
erate an ensemble of configurations, and a complex phase which is folded into the observables. 
This approach has two major drawbacks: the sign problem and the overlap problem. In order to 
address the overlap problem, a method which employs the canonical partition function has been 
proposed [|I|]. 

To generate ensembles using this method, the projected determinant for a fixed net quark 
number (number of quarks minus the number of anti-quarks) needs to be calculated using the 
Fourier transform of the fermionic determinant [^, ||]. To evaluate it, a discrete version of the 
Fourier transform was used |^] which was shown to be accurate [J3|] for quark numbers as large 
at 15. However, when simulating larger lattices, to achieve the same densities we need to use 
a proportionally larger number of quarks. As it turns out, the discrete Fourier transform becomes 
unstable when we use quark numbers larger than 20. To overcome this problem, we have developed 
a winding number expansion method (WNEM), a method based on the expansion in terms of quark 
loops. Using this approach, we can evaluate reliably the projected determinant for large quark 
numbers. 

2. Canonical approach 

The canonical partition function for lattice QCD is related to the grand-canonical partition 
function via the fugacity expansion 



where k is the net number of quarks and Z c is the canonical partition function. Thus, the canonical 
partition function can be written as a Fourier transform of the grand-canonical partition function 



In this paper, we will consider only the case of two degenerate flavors. If we denote M as the quark 
matrix for one quark flavor, we have 



Z(V, 



T^) = ^Z c (V,T,k)e^ T 



(2.1) 



k 




(2.2) 




(2.3) 



/ 



Redet k M 2 [U}\ det k M 2 [U] 
detM 2 [U] | Re det k M 2 [U] \ 



where 




(2.4) 



is the projected determinant with a fixed net quark number k. 
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3. Discrete Fourier transform instability 



The simplest way to compute the Fourier transform, Eq. ( |2.4| ), is to replace the continuous 
transform with its discrete version, 



det k M 2 [U] 



i N-l 



ik<pi 



N 



detM 2 [U, ty] 



<t>j 



7=0 



2nj 
N ' 



(3.1) 



Unfortunately, the Fourier transform becomes unstable when the quark number k is large as we can 
see from Table [j] This is a well known problem; it is a consequence of the fact that the higher 
Fourier components are the result of very delicate cancellations in a a sum of alternating terms. 



k 


3 


6 


9 


12 


15 


N=51 


2212.21 


247.601 


-22.8783 


-4.53755 


-0.233997 


N=102 


2212.21 


247.601 


-22.8783 


-4.53755 


-0.233997 


N=204 


2212.21 


247.601 


-22.8783 


-4.53755 


-0.233997 


k 


18 


21 


24 


27 


30 


N=51 


-0.00545724 


-0.0000602919 


6.70879E-7 


6.70879E-7 


-0.0000602919 


N=102 


-0.005458 


-0.0000631063 


-8.98294E-7 


1.56917E-6 


2.81435E-6 


N=204 


-0.005458 


-0.0000634881 


-1.66312E-6 


-2.83726E-7 


6.42123E-6 



Table 1: The projected determinant for a particular configuration with different choices of N. For k larger 
than 20 the results differ for different choices of N, signaling a numerical instability. 

For our studies using 6 3 x 4 lattices, to investigate the phase diagram at non-zero baryon 
densities, we need quark numbers of the order k ~ 30 which corresponds to 10 baryons. To evaluate 
the projection accurately we need a better method. 



4. Winding number expansion method 

The basic idea of the new method is to use the Fourier transform of logdetM(£/, 0) instead of 
the Fourier components of detM(£/, (p). Using an approximation based on the first few components 
of logdetM(£/, 0) we can then analytically compute the projected determinant. The success of the 
method can be traced to the fact that the Fourier components of logdetM(£/,0) are exponentially 
smaller with increasing order in the expansion 1 . This is why we can approximate the exponent 
very accurately with few terms which, in turn, allows us to evaluate the Fourier components of the 
determinant precisely. 

To see how this works we look at the hopping expansion of the logdetM(t/, </>). We start by 
writing the determinant in terms of its log 

detM(t/,0) =exp(TrlogM(t/,0)). (4.1) 

It is well known that TrlogM corresponds to a sum of connected quark loops. We separate all 
these loops in classes in terms of the net number of times they wrap around the lattice in the "time" 

'This property is also noted in a similar context in [n]. 
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Figure 1: Schematic of winding number expansion on lattice. 

direction (see Fig. [[]). We have then: 

TrlogM(I/,0) = £ L(U^) 

loops 



(4.2) 



= A (C/) + \Ee**W n (U) + e-^W,? (£/)], 
where n is winding number of quark loops and W n is the weight from the contribution of all the 



quark loops which have fixed winding number n in time direction. Eq. ( |4.2[ ) can be re-written as 
TrlogM(t/,0) =A (t/) + [£^W, ! (t/) + ^W„ t ([/)] (4.3) 

n 

= Aq (U) + Y, A n cos + 5„) , 



where A„ = 2|W„ | and <5„ = arg W„ are independent of 0. Using Eq. ( fkl| ) and Eq. ( p3| ) we get 

detM(f/,0) = e A o+Aicos(^+5i)+A 2 cos(20+5 2 )+ _ (44) 

To first order in the winding number expansion we have 

detAf(I/,0)„ =1 =<A+Ai cos^+SO _ (4.5) 
The Fourier transform can now be computed analytically; we have 



A0+A1 cos(0+5i) _ Ao+ikSj 



10 27t 

where is Bessel function of the first kind. 



(4.6) 
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For higher orders in the winding number expansion, we compute the Fourier transform using 
the Taylor expansion: 

f ^ ^ e -i^ g Ao+Aicos(0+5i) g A 2 cos(20+&)+A3Cos(3^+5 J )+... 

Jo 2n 

= I ^^M o+AiCOS ^ +5i »(1+A 2 cos(20 + 52) + 1a|cos(20 + 52) 2 + ...) x 
JO 27T 2 ! 

(1 +A 3 cos(30 + S3) + jfi cos(30 + S3) 2 + ...) x ... 

= Coolk(M)+C + o\Ik+l(Ax)+C-oi\-\{M)+C+02h+2{A\) +C_02l*-2(Al) + ... (4.7) 

In leading to the final result in Eq. (|i~7|), we have used the property that cos(«0 + 8 n ) = 
[ e '{n<l>+8n) _)_ e -'(' 5 0+ 5 «)] / /2 and the Fourier transform with two exponentials, according to Eq. (|4~6|), 
gives the Bessel functions 4± n (Ai). The coefficients coo,c + o n ,c_on are the sum of terms for the 
Bessel functions h{M),h=\ (Ai),4_„(Ai). 

Using Eq. ( p~7| ) and the recursion relation for the Bessel function, (A) = (A) (A), 
the winding number expansion can be extended easily to higher orders. 



5. Numerical tests 



As we mentioned before, the success of the method rests on the assumption that we can get 
a very good approximation of the exponent using only a few terms in the Fourier expansion. In 
Fig. H), in the top row we plot TrlogM(0) evaluated at 204 different values of (j) in the interval 
[0,2tt] (the three plots are the same); the second row shows an approximation using 1, 3 and 6 
terms respectively. It is easy to see that the approximations are all very good - this is due to 
the fact that the first term in the expansion is much larger than the subsequent terms. To see the 
contributions of the higher terms in the third row, we plot the difference between the exact value 
(top row) and the approximation (second row). It is easy to see that the error of the approximation 
decreases very rapidly with the number of terms. Note also that the error is well described by a 
cosine function: this is due to the fact that the Fourier coefficients A„ decrease exponentially with 
increasing n. The error will be dominated by the first term that is not included in the approximation: 
A2 for the first column, A4 for the second and A7 for the third. 

For our tests we decided to keep a winding number expansion to the sixth order. From Fig. || it 
is easy to see that the error introduced by this approximation is of the order of 10~ 7 which is precise 
enough for our purposes. It is important to note that to determine the coefficients A n exactly we 
would need to evaluate logdetM(£/, 0) for all possible values of (p. However, we can approximate 
their values using a discrete Fourier transform. To determine the right number of points Nd for the 
discrete 0,- needed, we compute these coefficients using an increasing number of Nd- In Table ^| we 
present the coefficients determined for a particular configuration; it is easy to see that Nd = 16 is 
good enough to determine the coefficients with sufficient precision for our purpose. 

The final test, and the most important, is to determine how accurate this approximation is when 
computing the projected determinant. In Fig. || we compare the result of this approximation with a 
high order approximation that uses the simple discrete Fourier transform. The expectation is that as 
we increase the number of quarks the value of the projected determinant decreases exponentially. It 
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Figure 2: A comparison between the exact value for logdetM(£/, (j)) and its winding number expansion to 
various orders. The value of logdetM(f/, <j>) is shifted so that it averages to zero. 





Nd=16 


Nd=24 


Nd=36 


Nd=200 


AO 


2.357308E+02 


2.357308E+02 


2.357308E+02 


2.357308E+02 


Al 


-1.341400E+01 


-1.341400E+01 


-1.341400E+01 


-1.341400E+01 


A2 


2.820535E-02 


2.820535E-02 


2.820535E-02 


2.820535E-02 


A3 


4.135219E-04 


4.135043E-04 


4.134942E-04 


4.134755E-04 


A4 


2.148188E-04 


2.147950E-04 


2.147792E-04 


2.147547E-04 


A5 


2.641758E-05 


2.639153E-05 


2.637794E-05 


2.636227E-05 


A6 


2.289491E-06 


2.286772E-06 


2.291285E-06 


2.305249E-06 



Table 2: The expansion coefficients determined using a discrete Fourier transform. 



is easy to see that this is indeed the case for all approximations at low quark number, say for k < 6. 
However, as the quark number gets close to 30 the projected determinant calculated from the old 
approximation (the red points) flattens out - this is the onset of numerical instability. The new 
approximation method doesn't suffer from this problem and the projected determinant continues 
to decrease as the quark number is increased. We note that the winding number expansion method 
(WNEM) approximates well the results from the high precision discrete Fourier transform (Nd = 
204) in its stability range. We conclude that the WNEM approximation is stable and accurate for a 
much larger range of quark number projection. 

From Fig. |3] we also note that, as concluded above, to determine the coefficients A n with 
enough precision, we don't need to use a high order approximation. Both the curves WNEM(16) 
and WNEM(204) with Nj = 16 and 204 respectively are determined using a 6 th order WNEM. The 
only difference is the set of coefficients used: WNEM(16) uses the coefficients appearing in the 
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Figure 3: The projected determinant as a function of the quark number. We plot here the results of various 
approximation based on either the discrete Fourier transform or on the winding number expansion. 



first column in Table g, while WNEM(204) uses the coefficients in the last column. We find that 
the resultant projected determinants differ only by 10~ 6 in relative errors for the quark number k 
up to 50. This is very important since it allows us to speed up the calculation tremendously without 
losing precision. 

To compare the merits of the WNEM with the discrete Fourier transform, we compute the 
values of the projected determinant calculated from the discrete Fourier transform with = 16 
only. In Fig.|| these results are labeled Discrete(16). It is easy to see that this approximation is only 
valid up to k = 6. This is to be compared to WNEM(16) that takes the same computational time 
but can approximate the projected determinant for k = 40 and beyond. 

One final test was to compare the 6 th order WNEM with the 1 th order WNEM. The results are 
plotted in Fig. ^. The relative error between the two is less than 10~ 4 for k up to 50. We conclude 
that a 6 th order WNEM using Nd = 16 is precise enough for the simulations for the present lattice 
(6 3 x 4 at j8 = 5.2 and quark mass K = 0. 158). 

6. Conclusion 

In this paper we discussed a novel approximation method for the projected determinant. The 
original approximation that employs a discrete Fourier transform becomes numerically unstable 
at large quark numbers which limits the region of the phase space that can be studied using the 
canonical approach. We find that our method based on the winding number expansion is numeri- 
cally stable and much faster which allows us to study higher baryon numbers. 

The winding number expansion method has been used successfully to carry out finite density 
simulations for QCD with two and four degenerate quark flavors [Q]. Our numerical checks show 
that for all these ensembles the approximation is accurate. 

Finally, we should mention that the conclusions we have drawn in this paper are valid for the 
fairly heavy quark mass we used (m % ~ lGeV). As the quark mass is lowered, we expect that more 
terms of A n will be needed in the expansion. 
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